Linear Scaling Solution of the Coulomb problem using wavelets 
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The Coulomb problem for continuous charge distributions is a central problem in physics. Powerful 
methods, that scale linearly with system size and that allow us to use different resolutions in different 
regions of space are therefore highly desirable. Using wavelet based Multi Resolution Analysis we 
derive for the first time a method which has these properties. The power and accuracy of the 
method is illustrated by applying it to the calculation of of the electrostatic potential of a full 
three-dimensional all-electron Uranium dimer. 
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The theory of waveletsflU is one of the most impor- 
tant recent development in mathematics. It allows one 
to apply a multi-scale analysis to problems that exhibit 
widely varying length scalesa. Problems with this fea- 
ture abound in all fields of physics. The problem we 
want to address here is the classical Coulomb problem 
for a continuous charge distribution p, i.e. we want to 
solve Poisson's equation 



(1) 



under the constraint that the potential V vanishes at in- 
finity. This basic equation can be found in nearly any 
field of physics and it is therefore essential to have effi- 
cient solution methods for it. There are two important 
requirements for an algorithm that solves this problem. 
First it should scale linearly with the size of the charge 
distribution. Since in numerical applications the charge 
distribution is given on a grid a measure for the size of the 
charge distribution is the number of grid points necessary 
to represent it. This linear scaling property is of utmost 
importance because in many applications one needs grids 
consisting of a very large number of grid points. Second, 
it should allow for grids that are nonuniform, i.e. that 
have higher resolution in regions where this is required. 
For discrete charge distributions several algorithmsQ with 
these two properties exist and have become a standard for 
simulations of large coulombic and gravitational particle 
systems. For continuous charge distributions proposals 
have been put forward to map the continuous problem 
onto a discrete and|-jase the above mentioned algorithms 
for discrete systemsB. To the best of our knowledge, there 
exists however no linear scaling algorithm for nonuniform 
grids that can directly be applied to continuous charge 
distributions. If one constrains oneself to uniform grids 
and periodic boundary conditiojis, there are of course the 
well known Fourier techniquea3, that show a nearly lin- 
ear Nlog2{N) scaling with respect to the number of grid 
points N . Non-periodic boundary conditions can be im- 
plemented in the context of Fourier techniques only by 
cutting off the long range Coulomb potentialQ. Finite 



element methods allow nonuniform grids, but grid gen- 
eration and preconditioning pose severe problems. Using 
a basis of wavelet functions, we will present in this pa- 
per a method that scales strictly linear and allows for 
nonuniform grids. 

There are many families of wavelets and one has to 
choose the most appropriate one for a specific applica- 
tion. A widely used family are thcrcompactly supported 
orthogonal wavelets of DaubechiesH. The orthogonality 
property is convenient if one has to expand an arbitrary 
function in a basis of wavelets. Their disadvantage is 
that they are not very smooth, i.e. only a small num- 
ber of derivatives is continuous. One can however con- 
struct families of nonorthogonal wavelets that are much 
smoother. In general the mapping from the numeri- 
cal values on a grid to the expansion coefficients of the 
wavelet basis is rather complicated and slow for biorthog- 
onal wavelet^. An exception are the second generation 
interpolating waveletsQ, whose special properties allow us 
to do this mapping easily. In addition they give rise to a 
particularly fast wavelet transfornJS. 

Wavelets have already been successfully applied in sev- 
eral areas of physical. In the context of electronic struc- 
ture calculations a fairly small wavelet basis can describe 
the widely varying scales of both core and valence elec- 
tronalj. Self-consistent electronic structure calculations 
have also been donetil. In these self-consistent calcula- 
tions the solution of Poisson's equation was however done 
by traditional Fourier techniques. 

Let us now, briefly review the theory behind biorthogo- 
nal waveletsH. As in ordinary orthogonal wavelet theory 
there are two fundamental functions, the scaling function 
(j> and the wavelet ^jj. In the biorthogonal case there are 
however still the complementary scaling function cf) and 
the complementary wavelet ip. Each scaling function and 
wavelet belongs to a hierarchical level of resolution. By 
analysing a function with respect to these different lev- 
els of resolution one can do a so-called Multi Resolution 
Analysis (MRA). The space belonging to a certain level 
of resolution k is spanned by all the integer translations 
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of the scaling function (f>i^k{x) oc (t){{^)'^x — i). Any func- 
tion can be expanded within this level of resolution. 

f{x) « ^Si,fc(/)i,fc(a;) , (2) 

i 

Since the scaling functions and their complementary 
counterparts are orthogonal < 4>k,i(x) \ 4>k,jix) >= 5ij, 
the expansion coefficients Si^k are given by 

Si,k =< f{x) I 4)i,k{x) > . (3) 

The expansion |^ becomes more accurate if one goes to 
a higher level of resolution, i.e. if one decreases k and it 
becomes exact in the limit A: — > — cx). This is an impor- 
tant feature because it allows us to improve systemati- 
cally the numerical accuracy in very much the same way 
as it is done with a basis of plane waves. In numerical 
application it is of course not possible to take this limit, 
and we will therefore denote the finest level of resolution 
that is used in the calculation by A: = 

The scaling function satisfies a refinement relation 

(f>j,k{x) = '^hi-2j(l)i,k~i{x) (4) 

i.e. each scaling function of a lower resolution level can 
be expressed as a linear combination of higher resolution 
scaling functions. It is obviously not possible to express 
a scaling function of higher resolution by a linear com- 
bination of lower resolution scaling functions only. One 
can, however, write down such an expression if one still 
includes the wavelets tpi^k{x) oc ^{{^)'^x — i) 

(j)j^k-i{x) = ^hi^2j(f>Lk{x) + ^gi^2j'ipi,k{x) (5) 

The wavelets at level k thus reintroduces the resolution 
that is lost as one goes from level k to level fc -I- 1 scaling 
functions. These transformation properties among the 
scaling functions and wavelets give rise to the wavelet 
transform. The wavelet expansion coefficients di^k are 
then either defined by this transform or equivalently by 
an expression analogous to Eq. ^. 

Si,k+1 = ^ /ji-2jSj,fe , di^k+1 — '^^gi-2jSj,k (6) 
j 3 

j 

Eq. ^ is called a forward fast wavelet transform, Eq. 
being its inverse counterpart. If one has periodic bound- 
ary conditions for the data Si^k the wavelet transform 
is one-to-one transformation between Si^k and its spec- 
tral decompositions Si^k+i and di^k+i- To obtain a full 
wavelet spectral analysis the forward transform is ap- 
plied recursively. In consecutive steps the output data 



Si.k+i of the previous forward transform are the input 
data for the next transform. The size of the data set to 
be transformed is thus cut into half in each step. The 
total operation count is then given by a geometric series 
and scales therefore strictly linear. A full wavelet syn- 
thesis consists in the same way of a sequence of inverse 
transforms and gives back the original data set Si^k- The 
coefficients hi and gi and their complementary counter- 
parts hi , gi are filters of finite length 2m and can be 
derived from the MRA requirementsu. -The 8-th order 
lifted Lazy scaling function and waveletO that were used 
in this work are shown in Fig.l. Because it can repre- 
sent polynomials up to degree 8 exactly, the expansion 
coefficients with respect to the wavelets di^k decay very 
rapidly for any smooth function with decreasing k. 
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FIG. 1. The scaling function (full line) and wavelet (dashed 
line) used in this work 

To do a multidimensional MRA, we use a scheme de- 
scribed by Daubechiesa. Even though all this work was 
done in the three-dimensional case, we will illustrate the 
principle just for the two-dimensional case. The space of 
all scaling functions of resolution level k is given by 

<l>i,3,k {x,y) ^ (l>iM {x)(t)j,k(y) (8) 

The wavelet space is again defined as the space that rec- 
ompensates for the resolution that is lost by going up 
one level in the scaling functions space. Using Eq. || one 
obtains three kind of terms for the wavelet space 

i'°'j,k{^,y) = M^)^jMy) (9) 
^llk(.x,y) = ^^Ax)(bJAy) (10) 
i^ljA^'y) = '^i,k{x)ipj,k{y) (11) 

Let us now explain how to solve Poisson's equation in 
wavelets. Expanding both the charge density and the 
potential in Eq |l| into scaling functions 

p{x, y,z) pi^j (f)i^k (x) (f>j,k iy) (12) 

V{x,y,z) = ^ V,j (j)i^k (x) (jijM iy) (13) 

one obtains the following system of equations. 

X] L^,j;t^,'^ "^M,^ = Phj (14) 
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where 



(15) 



matrix L'' 

Tk 

^il.i2;jl.j2 



Since the scahng functions have a finite support, the 
is a sparse matrix and its noiizero elements 
can be calculated analyticalljO. The natural 
boundary conditions for this scheme are periodic bound- 
ary conditions. As we stressed in the introduction, we 
however want to solve Poisson's equation ^ with non- 
periodic boundary conditions. As is well known, bound- 
ary aff'ects vanish whenever the boundary is sufficiently 
far away. Thus, one could in principle obtain natural 
boundary conditions (i.e. V{r) ^ if r — > oo) within 
arbitrary precision if one uses a sufficiently large peri- 
odic box. Since the electrostatic potential decays fairly 
slowly a very large box is required and the numerical ef- 
fort would be tremendous if one uses equally spaced grids 
within this huge periodic computational box. Far away 
from the charge distribution the variation of the potential 
is however small and less resolution is needed. The key 
idea is therefore to use a set of hierarchical grids as shown 
in Fig. 2. where the resolution decreases as on goes out of 
the center. Expressed in the terms of wavelet theory this 
means that on the highest (periodic) level we have a ba- 
sis of scaling functions. Resolution is then increased by 
adding wavelet basis functions near the center. By doing 
this repeatedly we obtain increasing resolution towards 
the center as shown in Fig. 2. 

Up to now the motivation for introducing grids of dif- 
ferent resolution was to handle the natural boundary 
conditions. Additional levels of resolution can however 
be introduced to handle charge distributions that have 
different length scales and require therefore higher res- 
olution in some parts of space. The theory of wavelets 
gives us also enough flexibility to increase the resolution 
not only around one center but around any number of 
centers in the computational box. 
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FIG. 2. A hierarchical multi resolution grid of the type used 
in this work. For simplicity only three levels of resolution are 
shown 



A full wavelet synthesis step can be done straightfor- 
wardly in this hierarchical grid setting. Any wavelet can 
be decomposed into scaling functions and therefore one 
can calculate the scaling function coefficients at any level 
of resolution and for any point in the computational vol- 
ume. If one calculates these scaling functions for high 



resolution levels in a region of low resolution, one ob- 
tains however a highly redundant data set. To do a full 
wavelet analysis that brings back the original spectral 
decomposition data, it turns out that one needs actually 
a slightly redundant data set. In order to calculate the 
wavelet coefficients for a wavelet at a boundary to a lower 
resolution region, one needs the scaling function values 
corresponding to this higher resolution also in a strip of 
width m in the lower resolution region. A schematic dia- 
gram of a full hierarchical wavelet analysis and synthesis 
is shown in Fig. 3. 
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FIG. 3. A schematic representation of a multi hierarchy 
wavelet synthesis. Data regions denoted by s,dP^ ,(1^" and d^^ 
contain expansion coefficients for basis functions of the type 
given by Eq. ^ ^ ^ and |l^ respectively. CP stand for a copy 
step where one puts an additional layer of zeroes around the 
data set. SYN denotes a one level wavelet synthesis step. One 
starts the process at the coarsest (periodic) level and proceeds 
down to the finest resolution level. To do a multi hierarchy 
wavelet analysis one proceeds back up reversing all the copy 
operations and replacing the single level synthesis steps by 
analysis steps. 



In this mixed representation, where one has scaling 
functions at the highest periodic level and wavelets all 
the refinement levels, the structure of the Laplace oper- 
ator is much more complicated since one has coupling 
between all the hierarchical levels. An elegant way to 
cope with this additional complexity is the so-called non- 
standard oppcator form proposed by Beylkin, Coifman 
and Rokhlii£3, which allows us to incorporate this cou- 
pling by a sequence of wavelet transforms (Fig. 3), that 
are interleaved with the application of a simple one-level 
Laplace operator. For this one-level Laplace operator 
only the matrix elements of the Laplace operator among 
scaling functions and wavelets on the same resolution 
level, but not between different levels of resolution are 
needed. 
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Mathematically the nonstandard operator form is a 
telescopic expansion of the Laplace operator in the scal- 
ing function basis at the finest level LF* . If we define 
projection operators Pk and Qk, that project the whole 
space into the space of scaling functions and wavelets at 
the k-th level as well as their complementary counter- 
parts Pk and Qfe, they satisfy 



Pk = Pi 



k+l 



Qk+1 



Pk = P« 



k+l 



!k+l 



(16) 



and we may write 



L'^ = Pk L" Pk - {Qk+i + Pk+i) {Qk+i + Pk+i) 



= L 



k+l 
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L 



k+l 
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k+l 



(17) 



where L'^j^, ^s'I) i ^^ds Laplace operators at the 
{k + l)th level representing the coupling of wavelets 
with wavelets, wavelets with scaling functions and scal- 
ing functions with wavelets. Applying Eq. ^ recursively 
for k = 0,1,..., one obtains the nonstandard operator 
form. 

In the basis of the wavelet functions at different reso- 
lution levels a simple diagonal preconditioning scheme is 
very efficient and we were able to reduce the residue by 
one order of magnitude with only 3 iterations. 

To demonstrate the power of this method we applied 
it to a problem that can hardly be solved by any other 
methods, namely the potential arising from the nucle- 
onic and electronic charge distribution of a fully three- 
dimensional all-electron Uranium dimer. The charge 
distribution of the nucleus was represented by a Gaus- 
sian charge distribution with an extension of — ^ atomic 
units. Since the valence electrons have an extension 
which is of the order of one atomic unit, we have length 
scales that differ by more than 3 orders of magnitude. 
As can be seen from Fig. 4, the potential also varies by 
many orders of magnitude. Using 22 hierarchical levels in 
our algorithm, we can represent resolutions that differ by 
7 orders of magnitude and we are able to calculate the 
potential with at least 6 significant digits in the whole 
region from the nucleus to the valence region. In order 
to be able to determine the error we actually first fitted 
the electronic charge distribution by a small number of 
Gaussians, whose exact potential can be calculated ana- 
lytically. This rather crude charge density was then used 
in all the calculations. 

We also applied the method to clusters containing 
several CO molecules that were described by pseudo- 
potentials. In this case there is only one length scale 
associated with the charge distribution and it is possible 
to reduce the number of grid points on higher levels such 
that the total amount of work increases only slightly with 
additional hierarchies. We were able to calculate the po- 
tential corresponding to the non-periodic boundary con- 
ditions with 8 significant digits. 

We thank Mike Teter and Leonid Zaslavsky for bring- 
ing the beauty and usefulness of wavelets to our atten- 
tion. Jurg Hutter pointed out several essential references 



on wavelets to us. We acknowledge the interest of O. K. 
Andersen, O. Gunnarson, and M. Parrinello. 
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FIG. 4. The potential V (full line) and the numerical error 
AV^ (dashed line) for an Uranium dimer as a function of the 
distance from of right hand side nucleus. The left panel shows 
both quantities in the direction of the left nucleus, the right 
panel shows them in the opposite direction. Both distances 
are given on a logarithmic scale. 
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